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Abstract: 

In this contribution a broad overview of the methodologies of cosmological iV-body simulations and a 
short introduction explaining the general idea behind such simulations is presented. After explaining 
how to set up the initial conditions using a set of N particles two (diverse) techniques are presented 
for evolving these particles forward in time under the influence of their self-gravity. One technique 
(tree codes) is solely based upon a sophistication of the direct particle-particle summation whereas the 
other method relies on the continuous (de-)construction of arbitrarily shaped grids and is realized in 
adaptive mesh refinement codes. 

Keywords: cosmology: theory - cosmology: dark matter - methods: n-body simulations - methods: 
numerical 
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1 Introduction 

The purpose of cosmological simulations is to model 
the growth of structures in the Universe. They have a 
long history and numerous applications. These simula- 
tions play a very significant role in cosmology because 
they can be considered as an "experiment" to verify 
theories of the origin and evolution of the Universe. 

The Universe is believed to have started with a Big 
Bang in which - or more precisely: shortly after which 
- tiny fluctuations (in an otherwise homogeneous and 
isotropic space) were imprinted into the radiation and 
matter density field. To understand how the Universe 
evolved from that early stage into what we observe 
today (i.e. stars, galaxies, galaxy clusters, ...) one 
needs to follow the evolution of those density fields us- 
ing numerical methods as soon as they turn non-linear. 
Therefore, the approach to cosmological simulations is 
actually twofold: firstly, one needs to generate the ini- 
tial conditions according to the cosmological structure 
formation model to be investigated and secondly, the 
initial density field (sampled by the particles) needs to 
be evolved forward in time using an iV-body code. 

In all such codes the evolution is simulated by fol- 
lowing the trajectories of particles under their mutual 
gravity. These particles are supposed to sample the 
matter density field as accurately as possible and a 
cosmological simulation is nothing more (and noth- 
ing less) than a simple and effective tool for investi- 
gating non-linear gravitational evolution. There are 
two constraints on a cosmological simulation though: 
a) the correct initial conditions and b) the observa- 
tion of galaxies, galaxy clusters, large-scale structure, 
voids, etc. Simulations are hence trying to bridge the 
gap between observations of the early Universe (i.e. 
anisotropies in the Cosmic Microwave Background ob- 
served as early as 300000 years after the Big Bang) 
and the Universe as we see it today. 

The first application of iV-body methods in astro- 



physics was in the simulations of star clusters using as 
little as a handful of particles (Aarseth 1970). During 
the 1970's more simulations of galaxy clustering were 
performed using what is called PP (particle-particle) 
methods (Peebles 1970) and not earlier than 1981 the 
first cosmological simulations using more than 20000 
particles became feasible (Efstathiou et al. 1981). 

Until now the methods have been continuously re- 
fined to allow for more and more particles while simul- 
taneously resolving finer and finer structures. Today 
it is standard to run a cosmological simulation with 
millions of particles in a couple of days on large super- 
computers or even clusters of PC's (cf. Gill, Knebe & 
Gibson 2004) . These simulations can resolve the orbits 
of satellite galaxies within dark matter haloes spanning 
five orders of magnitude in mass and a spatial dynam- 
ical range well above 30000. 

In this contribution I would like to focus on two nu- 
merical techniques in particular presenting their method- 
ologies, advantages and shortcomings when being com- 
pared with each other. I need to stress though that 
all cosmological simulations are based on the assump- 
tion that the Universe mainly consists of dark matter 
interacting merely gravitationally. Baryonic matter, 
which only accounts for about 15% of the total mass, 
is accounted for in such simulations only via its grav- 
itational effects, too. Even though todays simulation 
methods and computer technology have become suf- 
ficiently sophisticated as to allow for hydrodynamical 
processes to be included, the particulars of such im- 
plementations, however, are beyond the scope of this 
contribution. 

Fig. Q depicts the conceptual ideas behind (cos- 
mological) iV-body simulations: starting from initial 
seed inhomogeneities superimposed onto a homoge- 
neous and isotropic background the matter field is evolved 
forward in time. This evolution depends on the cos- 
mological model under investigation and is performed 
using an iV-body code. Snapshots of the simulation 
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at various times are recorded and then analysed and 
compared to observational data to verify and falsify 
theories of structure formation and evolution. 

2 The Initial Conditions 

The most common way to set up initial conditions for a 
cosmological simulation is to make use of the Zeldovich 
approximation to move particles from a Lagrangian 
point <f to a Eulerian point x (e.g. Efstathiou et al. 1985) 

X = q-D(t)S(q) , (1) 
where D(t) describes the growing mode of linear fluctu- 
ations and S(q) is the 'displacement field'. The initial 
Lagrangian coordinates q are usually chosen to form a 
regular, three-dimensional lattice. 

The displacement field S(q) is related to a pre- 
calculated input power spectrum of density fluctua- 
tions, P(k), which in turn depends on the cosmological 
model under consideration 

§(q)=V q $(q), $(g) = ^ a s cos(fc-^+fej sin(fc-g) , 

k 

(2) 

where the Fourier coefficients and 6j are linked to 
P(k) and are given as 

H = R^VW), h = R^Vm- (3) 

Ri , R2 are (Gaussian) random numbers with zero mean 
and dispersion unity. 

An example for the Zeldovich approximation "at 
work" can be found in Fig. |2Jwhere a slice through the 
initial conditions for a standard ACDM simulations at 
redshift z = 50 is shown. This figure nicely demon- 
strate the decomposition of the density perturbations 
as waves. 

3 Solving Poisson's Equation 

When used to model the dynamics of a collisionless 
system such as dark matter, an TV-body code aims 
at simultaneously solving the collisionless Boltzmann 
equation (CBE) 

df ( df 9$ df A 

i=l 

and Poisson's equation 

V 2 $(r) = A-nGp{r) . (5) 

The CBE Eq. Q is solved by the method of character- 
istics (e.g., Leeuwin, Combes & Binney 1993). Since 
the CBE states that / is constant along any trajectory 
{f(t),v(t)}, the trajectories obtained by time integra- 
tion of TV points {ri, Vi} sampled from the distribution 
function / at time t = tmitiai form a representative 
sample of / at each time t. 



Hence the problem reduces to solving Poisson's 
Eq. (S) for a set of TV particles and advancing them 
forward in time according to the equations of motion 
derived from the system's Hamiltonian H (remember 
that Eq. 0} can be written as df/dt + [f,H] = 0). 
The details of the time-integration of the equations of 
motions are going to be explained later on in Section[I] 
though; in this Section I like to focus on the gravity 
solver. 

Currently there are two commonly used approaches 
for deriving the potential from Poisson's equation: a) 
tree codes rely on a direct particle-particle summation, 
and b) PM (particle-mesh) codes utilize a numerical 
integration of Eq. @ on a grid. 

3.1 Tree Codes 

3.1.1 The Pre-Requisites 

The particle-particle (PP) method upon which tree 
codes are based assumes that the particles are ^-functions 
and hence the density field (rhs of Poisson's equa- 
tion Eq. JSJ) reads as follows: 

JV 

p{r) = y^mi5(f- ft) , (6) 

1=1 

where TV is the total number of particles in use. 

3.1.2 The Forces 

Combining Eq. © with Eq. the analytical solution 
for the force F at particle position is given by: 

3^ 

But as we are interested in deriving the force at every 
single particle position, the PP method scales like TV 2 
(TV summations, each over (TV — 1) particles). There- 
fore, a (straightforward) PP summation appears not 
to be feasible for evolving a set of TV particles under 
their mutual gravity, not even on the largest super- 
computers available nowadays! One needs to bypass 
the increase in computational time for large numbers 
of particles with a more sophisticated treatment when 
calculating the forces. One way of achieving this is 
to organize the particles in a tree-like structure: par- 
ticles located "far away" from the actual particle (at 
which position we intend to calculate the force) can 
be lumped together as a single - but more massive - 
particle. This tunes down the number of calculations 
dramatically. 

The idea of a tree code is sketched in Fig. El The 
particles are organized in a tree-like structure based 
upon a cubical decomposition of the computational do- 
main. Consequentially, for each particle we "walk the 
tree" and add the forces from branchings that need no 
further unfolding into finer branches according some 
pre-selected "opening criterion" . 

One publicly available tree code is called GAD- 
GET 1 (Galaxies with Dark Matter and Gas intEracT) 

1 GADGET can be downloaded from at this web address 
http : / /www . mpa-garching . mpg . de/gadget 
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and I refer the reader to a more elaborate discussion 
of this technique to its reference paper by Springel, 
Yoshida & White (2001). 

3.1.3 Force Resolution 

In order to avoid the singularity for r» = f] in Eq. J7J 
one needs to set a limit on the minimal allowed spatial 
separation between two particles. This can be achieved 
by introducing a (fixed) scale, i.e. the softening param- 
eter e: 

This softening is closely related to the overall force 
resolution of the simulation and an elaborate discus- 
sion of it can be found in Dehnen (2001). 

3.2 Particle-Mesh Codes 

3.2.1 The Pre-Requisites 

Another way for obtaining the forces is to numerically 
integrate Poisson's equation ©. This method, how- 
ever, demands the introduction of a grid in order to 
define the density and hence the name particle-mesh 
(PM) method. The grid is usually of a regular (cubic) 
shape with L x L x L cells where each cell is identi- 
fied by the index triplet (i,j,k). The forces are then 
calculated according to the following scheme: 

1. assign all particles to the grid to get Pij,k 

2. solve Poisson's equation V 2 (f>i,j,k = iirGpi,i,k 

3. differentiate to get forces Fi,j,k = —^4>i,j,k 

4. interpolate Fi,j,t back to particle positions 

3.2.2 The Forces 

With this scheme most of the time is spent in step 2 
and the most common way to solve Poisson's equa- 
tion on a grid is to make use of FFT's (Fast-Fourier- 
Transforms). The analytical solution to Poisson's equa- 
tion is given by the integral 

$(r) = J G(r-f)p{r)dr (9) 

where G(x) = ~x/x 3 ^ 2 is the Green's function of Pois- 
son's equation. This integral can readily be evaluated 
in Fourier-space, i.e. 

$ = G p (10) 

where $, G, and p are the Fourier transforms of the 
respective variables. 

The PM approach proves to be exceptionally fast 
outperforming any tree code. 

There are of course other techniques than the use 
of FFT's available to numerically solve Poisson's equa- 
tion but the utilisation of FFT's is the most common 
approach as it appears to be the fastest. 



3.2.3 Force Resolution 

The most severe problem with the PM method is the 
lack of spatial resolution below two grid spacings. Whereas 
tree codes require the introduction of a softening length 
to avoid the force singularity for close encounters of 
particles PM codes suffer from the opposite problem. 
Gravity is an attractive force and hence the particles 
flow from low density regions into high density regions 
amplifying primordial density fluctuations. This leads 
to an excess of particles in certain cells whereas other 
cells are becoming more and more devoid of matter. 
But as the spacing of the grid introduces a (smooth- 
ing) scale particles closer than about two cell distances 
do not interact according to Eq. JJJJl anymore. While 
we had to introduce a force softening for tree codes 
to avoid two-body interactions and make the simula- 
tion collisionless, respectively, the PM method natu- 
rally admits such a smoothing. However, we are left 
with the situation where we can not resolve structure 
formation on scales of (and below) roughly the cell 
spacing of the grid! 

This is a major problem and the most obvious way 
to overcome it is to introduce finer grids in regions of 
high density. These grids though need to freely adapt 
to the actual particle distribution at all times and 
hence navigating such complex grids through computer 
memory is a very demanding task. One of the freely 
available adaptive mesh refinement (AMR) codes is 
MLAPM 2 (Knebe, Green & Binney 2001). 

The mode of operation of this AMR technique can 
be viewed in Fig. 0] where a slice through a standard 
ACDM simulation is presented. The left panel shows 
the distribution of particles whereas the right panel 
indicates the adaptive meshes used to obtain a solution 
to Poisson's equation. 

It needs to be stressed though that the use of ir- 
regularly shaped grids inhibits FFT's. Another tech- 
nique for solving Poisson's equation needs to be sought 
such as, for instance, multi-grid relaxation (cf. Knebe, 
Green & Binney 2001). 



3.3 Hybrid Methods 

There are, of course, various other techniques for si- 
multaneously being time efficient and having a credi- 
ble force resolution. One possibility is realized in the 
so-called P 3 M technique (i.e. Couchman 1991) where 
a combination of PP and PM provides the necessary 
balance between accuracy and efficiency: the force as 
given by the plain PM calculation is augmented by a 
direct summation over all neighboring particles within 
the surrounding cells. This gives accurate forces down 
to the scale provided by the softening parameter e 
again. Other examples, for instance, are Tree-PM 
(Bode & Ostriker 2003) and moving mesh (Gnedin 
1995) codes, but the details are well beyond the scope 
of this contribution. 



2 MLAPM can be downloaded from this web address 
http : //www . aip . de/People/AKnebe/MLAPM 
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3.4 Mass Resolution 

It still needs to be mentioned that a cosmological sim- 
ulation in practice only simulates a certain fraction 
of the Universe. This is what people refer to as the 
simulation box. However, to account for the fact the 
Universe is actually infinite one uses periodic bound- 
ary conditions: particles leaving the box on one side 
immediately enter the box again on the other side. 

Moreover, the size of the box also defines the mass 
resolution of the simulation. We are only using a cer- 
tain number of particles within a fixed region of the 
Universe. And as the density of the Universe is deter- 
mined by the cosmological model under investigation, 
each individual particle has a certain mass. This mass 
determines the mass resolution of that specific simula- 
tion. For instance, if we model the evolution of about 2 
million particles in a box with side length 25/i _1 Mpc 
using the ACDM cosmology (Qo = 0.3), each particle 
wcig hs about 6 • 10 s h' 1 M . Therefore we will not be 
able to resolve dwarf galaxies in that particular cos- 
mological simulation (Md war f > 10 7 /i _1 Mq). 

3.5 Comparison 

It only appears natural to ask the question which method 
is superior and how they compare, respectively. 

There is no straight forward answer as both meth- 
ods have their (dis-)advantages. Tree codes are based 
upon the assumption that the Universe is filled with 
particles of a certain size related to the softening e (cf. 
Section f3.1.3|l . Adaptive mesh refinement codes use a 
smoothed density field (which in turn also introduces 
an effective particle size) to obtain the potential and 
hence the force field. In both cases it is important to 
bear in mind that particles are only to be understood 
as markers in phase-space and should not interact on 
a two-body basis, i.e. one always intends to integrate 
the collisionless Boltzmann equation There are 

several studies investigating two-body interactions in 
such simulations and it can be confirmed that they are 
more prominent in tree codes (Binney & Knebe 2002). 
But as long as one complies with certain constraints 
on the numerical parameters (cf. Power et al. 2003) 
such effects can be minimized. 

There are several studies comparing tree and AMR 
codes both in efficiency and accuracy (cf. Frenck et al. 1999, 
Knebe et al. 2001, O'Shea et al. 2003) but in the end 
it all comes down to a "question of taste" . Both tech- 
niques are well enough developed to successfully model 
the formation and evolution of cosmic structures. 



4 Newtonian Mechanics in an 
Expanding Universe 

Even though solving Poisson's equation is the heart 
and soul of every iV-body code, it is also important 
to accurately update particle positions and velocities, 
i.e. integrating the equations of motion. And as the 
Universe is expanding it is convenient to introduce co- 
moving coordinates: 



x = r/a(t) (11) 

where a(t) is the cosmic expansion factor. 
The (comoving) Lagrangian is given by 



r l 2.2 $ 

L = ~a x , 

1 a 

which leads to the canonical momentum 



Hamilton's equations are therefore 



Ax 

"dT 



„2 



(12) 



(13) 



(14) 



dp V$ 

At a 
These equations can be discretized and integrated us- 
ing a second-order accurate scheme as follows: 



1+1/2 = X n +p„ 



t+At/2 



At 

„2 



Pn+i = p n - V$(a;, l+ i/2 



t+At 



At 



(15) 



X n + 1 = X n + 1 / 2 +Pn + l 



t+At 



At 



t+At/2 



where the integrals can be evaluated analytically as 
they depend only on the cosmology. 

This modified leap-frog scheme only needs to store 
one copy of the positions and velocities whereas other 
integrators as, for instance, Runge-Kutta consume more 
memory. 

5 Final Remarks 

A broad overview of the concepts behind cosmologi- 
cal iV-body simulations has been presented. After a 
short introduction explaining the needs for such sim- 
ulations I explained how to actually set up the initial 
conditions. There are numerous techniques for evolv- 
ing that set of N particles forward in time under the 
influence of its own self-gravity alone, but I focused on 
two (diverse) methods, namely tree codes and adap- 
tive mesh refinement codes. Their mode of operation 
has been introduced and their limitations pointed out 
and compared with each other. 

This contribution can only be understood as a very 
brief and general introduction to the ideas behind cos- 
mological iV-body simulations. It is far from being 
complete and exhaustive and I refer the reader to more 
elaborate review articles such as Bertschinger (1998) 
and the monograph "Computer Simulations using Par- 
ticles" by Hockney & Eastwood (1988). 
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Figure 1: Illustration of the idea driving A-body simulations. Initial conditions are being evolved forward 
in time modeling gravity alone. The outputs over time are then compared to observational data and 
the cosmological model adjusted accordingly. Image credit (mass map of C10024+164): European Space 
Agency, NASA and Jean-Paul Kneib (Observatoire Midi-Pyrenees, France/Caltech, USA) 
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Figure 2: Initial conditions according to the Zeldovich approximation (cf. Eq. Q). One clearly sees the 
decomposition of density perturbations as waves. 



Figure 3: Illustration of a tree code. The left panel shows the actual particle distribution and its cubical 
decomposition. The right panel is the tree corresponding to this distribution. 
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Figure 4: An example for adaptive mesh refinement. The left panel shows the particle distribution at 
redshift z = in a ACDM simulation. The right panel indicates the arbitrarily shaped grids invoked by 
the AMR code MLAPM to solve Poisson's equation. 



